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ABSTRACT 

We demonstrate that cosmic rays form filamentary structures in the precursors of 
supernova remnant shocks due to their self- generated magnetic fields. The cosmic-ray 
filamentation results in the growth of a long wavelength instability, and naturally cou- 
ples the rapid non-linear amplification on small scales to larger length scales. Hybrid 
magnetohydrodynamics-particle simulations are performed to confirm the effect. The 
resulting large scale magnetic field may facilitate the scattering of high energy cosmic 
rays as required to accelerate protons beyond the knee in the cosmic-ray spectrum at 
supernova remnant shocks. Filamentation far upstream of the shock may also assist 
in the escape of cosmic rays from the accelerator. 
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1 INTRODUCTION 

It is generally accepted that galactic cosmic rays are acceler- 
ated in supernova remnants. The process of diffusive shock 
acceleration (|Krymskii||1977| |Axford et al.||1977| |Bell||1978 



|Blandford fc Ostriker|1978j remains the most likely mecha- 
nism for producing and maintaining the observed spectrum 
There is now a growing wealth of observational evidence sup- 
porting this scenario. The detection of TeV gamma-ray emis- 
sion from nearby remnants confirms the presence of elec- 
trons, and possibly protons, with energies of at least 10 14 eV 
(e.g. |Hinton fc Hofmann|2009| ). In addition, high resolution 
observations of narrow non-thermal X-ray filaments at the 
outer shocks of several shell-type supernova remnants favour 
a model in which the highest energy electrons are produced 
directly at the shock, consistent with the predictions of dif- 
fusive shock acceleration (e.g. Vink & Laming 2003, Bamba 



et al.|2005| |Uchiyama et al.|2007[ ). These filaments also pro 
vide evidence for strong magnetic field amplification in the 
vicinity of the shock. The generation of strong magnetic tur- 
bulence is vital for the acc eleration of cosmic rays to the knee 
(~ 10 15 ' 5 eV) and above (Lagage & Cesarsky 



1983 



Bell & 



Lucek||200U 

While several mechanisms for amplifying magnetic 
fields to values in excess of the shock compressed interstel- 
lar fields have been proposed, those that result in the trans- 
fer of upstream cosmic-ray streaming energy to the mag- 
netic field are of greatest relevance for diffusive shock ac- 
celeration. The non-resonant mode first identified by |Bell| 
(20041 has been demonstrated to grow rapidly, however, 



the characteristic wavelength of the amplified field is pre- 
dominantly on a length scale shorter than the gyroradius 
of the driving particles. Under certain conditions, other 
short-wavelength instabilities may dominate (e.g 



Bret 



Riquelmc & Spitkovsky 2010; Lcmoinc & Pclletier 



2009 



2010 



Nakar et al.|2011 1. While such instabilities may be sufficient 



to explain the large magnetic field values inferred from ob- 
servations, in order to facilitate rapid acceleration to higher 
energies, it is necessary to generate field structure on scales 
comparable with the gyroradius of the highest energy parti- 
cles ( |Bell fc Lucek|2001|[Reville et al.|2008 l. 

The generation of large scale field structures has been 
investigated in the context of filaments, or beams, in |Bell| 
( |2005[ ), where a pre-existing profile for the cosmic ray dis- 
tribution was assumed. In this paper, we demonstrate that 
the distribution of relativistic particles is inherently non- 
uniform, and that filamentation occurs as a natural conse- 
quence of cosmic-ray streaming. A similar phenomenon oc- 
curs in laser plasmas whereby photon beams filament due 
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to thermal self-focusing in expanding cavities (e.g. Craxton 
|fc McCrory|1984| ). We show here that in the case of cosmic- 
rays, this process results in the growth of magnetic field on 
large scales. The development of the filamentation and large- 
scale field is investigated analytically in a two dimensional 
slab symmetric geometry, and verified using hybrid particle- 
MHD simulations. The non-linear feedback between ultra 
relativistic particles and the background plasma, and the 
resulting large-scale fields will have important implications 
for the acceleration of cosmic rays to energies above the knee 
in supernova remnants, and also their escape. 

The outline of the paper is as follows. In the next section 
we develop the analytic model that describes the cosmic-ray 
filamentation. It is demonstrated that this introduces a long 
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wavelength instability in the precursors of supernova rem- 
nant shocks. In section |3j we report on the numerical code 
used to investigate the instability, and present simulation 
results. The relevant time and length scales inferred from 
theory and observations are addressed in section [4] We con- 
clude with a discussion of the implications for cosmic ray ac- 
celeration and escape of cosmic rays upstream of the shock 
in supernova remnants in the context of filamentation. 



2 COSMIC RAY FILAMENTATION 

Within the diffusion approximation of shock acceleration 
theory, a first order anisotropy is introduced in the up- 
stream particle distribution as a result of the gradient in the 
isotropic part of the distribution. For a shock front propa- 
gating in the positive x direction, with velocity u s h, the re- 
sulting steady state test-particle solution is (e.g |Drury|1983[ ) 

F(x,p)=f (x,p) (l+3^- cos6>) (1) 

where fo(x,p) is the isotropic part of the spectrum as mea- 
sured in the upstream rest-frame, and 9 is the particle pitch 
angle with respect to the shock normal. It is generally under- 
stood that the current associated with the anisotropic part 
of the upstream particle distribution drives the growth of 
MHD instabilities. The resulting turbulent magnetic fields 
mediate the scattering that maintain the cosmic rays' quasi- 
isotropic distribution, ensuring a high probability that par- 
ticles repeatedly cross the shock before escaping. 

To investigate the behaviour in higher dimensions, it 
is convenient to use the Vlasov equation, which for ultra- 
relativistic particles can be written in the form 

df 
dp 

In a reference frame in which the upstream cosmic rays are 
isotropic, the distribution function f(x,p,t) is dependent 
only on the length of the momentum vector, not its direc- 
tion, and it is straightforward to show that the magnetic field 
term makes no contribution. It follows that the only force 
relevant for calculating the particle distribution is that due 
to the local electric field. On average, the upstream cosmic- 
ray distribution is isotropic in the rest frame of the shock 
1 McClements et al.||1996 |. Neglecting the bulk deceleration 
of the incoming plasma due to the cosmic-ray pressure gra- 
dient, the background plasma in this frame moves towards 
the shock with velocity u = — u s hX + 5u, where Su are the 
superimposed background fluid motions due to the cosmic- 
ray current, and x is the unit vector along the direction 
of the shock normal. Conservation of momentum dictates 
that these motions are small compared to the shock velocity 
\Su\ -C w s h, and to lowest order the local electric field is, in 
the ideal MHD limit, E — u s hX x B±. While this analysis is 
valid for all shock obliquities, we focus here on self-generated 
magnetic fields B± due to current-driven instabilities. 

To investigate the role of filamentation in the plane nor- 
mal to the direction of the cosmic-ray streaming, we neglect 
the cosmic-ray pressure gradient in the precursor, and re- 
strict the analysis to the case of slab symmetry in the x 
direction. Introducing the vector potential B = V x A, the 
local electric field is 



% + ■ V/ + e (E + v x B) 

at p 







(2) 



where the scalar potential is Ay = A ■ x. Inserting into Q, 
the distribution function, when observed in a reference frame 
in which the shock is at rest, evolves according to 



at p 



eV 1 



(4) 



where n s hA| | plays the role of th e effective electric field po- 
tential (e.g. Krall fc Trivelpiece||1973| section 8.17). On the 
slowly evolving MHD timescales, the cosmic-ray distribution 
will progress through equilibrium states 



= , where e = p — eu^A^ jc 



It follows that the phase-space distribution of the cosmic 
rays consists of surfaces of equal density on the momentum 
iso-surfaces e. Hence, if df /dp < 0, as is almost certainly the 
case, the cosmic-ray number density will be locally larger 
(smaller) in regions of positive (negative) Am. Specifically, 
if p 2> e« s h|A|i|/c, making a Taylor expansion, the number 
density as a function of position is 



r (y,z) = n + 



ettshAn 



87rp/ dp 



(5) 



where we have performed an integration by parts. Here fo is 
the unperturbed part of the spectrum, and no = J 4ivp 2 fodp 
the associated uniform number density. Note that the cor- 
relation with A11 is dependent on the choice of orientation. 
If the upstream instead was chosen to lie in the half plane 
x < 0, the density and vector potential would anti-correlate. 
In addition, the correlation is charge dependent, such that in 
the precursor, the electrons and protons will anti-correlate. 
Since the number density of non-thermal particles is very 
much less than that of the background plasma, on the length 
scales of interest, charge neutrality is always maintained. 

The growth of the magnetic field is driven by the result- 
ing cosmic-ray current. Transforming back to the upstream 
frame, from |5J, it follows that the cosmic-ray current is also 
a function of position 



3cr(y,z) =jo + x( A \\ ~ i A \\)) 



(6) 



where the additional term (A\\) has been added to conserve 
total particle number. Since lower energy particles are con- 
fined closer to the shock, the distribution is expected to fall 
off rapidly below a minimum momentum p m i n , where p m in(a;) 
increases with distance from the shock (e.g. |Eichler||1979| ). 
Assuming a power-law spectrum /o(p > Pmin) oc p~ 



2 2 
e u sh 



871-p/odp = 



2 2 
e n u ah 

PminC 



It has been implicitly assumed here that the distribution 
function is gyrotropic such that j CI x x — 0. This naturally 
holds on scales larger than the cosmic-ray gyroradius, but 
also on smaller scales provided small-angle scatterings on the 
background field are sufficiently frequent. All of these results 
have been verified using high resolution hybrid simulations 
(see section pjl. 

The role of a single cosmic-ray filament, or beam, with 
fixed cosmic-ray current has been previously investigated in 



Bell ( 2005 1. We have demonstrated here that the cosmic-ray 



E = u sh VA|i (y, z) 



(3) 



current is in fact filamentary in general. This will alter the 
growth of plasma instabilities, and in the next section we 
investigate this effect. 
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Figure 1. Illustration of the behaviour on different length scales. j CI is the cosmic ray current, and j r et = (V X B)//iq — j C i the 
return current carried by the background plasma. The small scale circles represent the expanding loops of magnetic field of a particular 
handedness. At early times (left), the magnetic field contains loops of both orientation having comparable strength, and the cosmic ray 
current is approximately uniform. Since only loops with favourable orientation can grow, in this example counter-clockwise, the cosmic 
rays are focused into these expanding loops, while clockwise loops contract. At later times (right) the large scale magnetic field Bgg on 
scale ro enclosing the smaller expanding loops, will have increased. 



2.1 Filament growth 

The amplification of magnetic fields in the precursors of su- 
pernova remnant shocks is a multi-scale problem. While the 
non-linear growth of magnetic field, driven by cosmic-ray 
currents has been demonstrated via numerical simulations 



(e.g. |Bell||2004| |Ohira et al.||2009| |Riquelme fc Spitkovsky 
|2009| |Stroman et al.||2009[ ), most simulations have focused 
exclusively on field amplification on small scales. This is pri- 
marily due to the fact that, on larger length scales, it be- 
comes necessary to include the dynamics of the cosmic rays 



( Lucek fc Bell|2000[ ) . We use the analysis of the previous sec- 
tion to self-consistently model the interaction between the 
background plasma and the cosmic rays. 



Following Bell ( 2005 1 , we analyse the MHD equations 



with an external cosmic-ray current, neglecting the role of 
pressure gradients and magnetic tension 

dB 



dt 
du 



V x (u x B) 

-jcr X B 



(7) 
(8) 



valid in the long wavelength approximation. Also, from Q, 
the vector potential satisfies 



OA 

~dt 



u x (V x A) 



(9) 



Returning to the 2D analysis of the previous section, 
i.e. zero gradient in the direction of cosmic-ray streaming, it 
follows that the parallel component of the vector potential 
is constant for a particular fluid element 

d^ll , , 

-dT = °- < 10 > 

This has a number of important consequences. From pjl, the 
cosmic ray current will increase in regions of large An. Con- 
sidering an idealised axisymmetric system, with maximum 
A\\ is centred on the origin, Be = —dA\\/dr > 0, at least 
locally. Hence, the resulting —j a x B force acts to push 
the plasma radially outwards. This spreads the region of 



large A\\ , thus focusing more cosmic rays into the filament, 
leading to a run-away instability. In a 2D slab symmetric 
geometry, this results in the spreading out of fiat table-top 
structures with large A\\ surrounded by regions of negative 
j4||, with large gradients in between. This is similar to the 



picture presented in Bell (2005) section 3, where the cosmic 



ray current was fixed, and the growth rate for the expansion 
of cavities was found to be 

rm= (h*) l/ ' , (ii) 

with r the radius of the cavity, and Bo the magnetic field 
strength on that scale. Here, the focussing of the cosmic rays 
into the cavities will enhance the growth rate as compared 
with the constant current case, since j CI is larger in the fil- 
aments. 

For growth on small scales the orientation of the mag- 
netic field must be favourable. Considering a field configu- 
ration such that at early times, an equal number of small- 
scale loops of both polarisations are randomly located within 
a circle of radius ro, as shown in Figure [l] the effect of 
the cosmic ray current is to expand loops of one orienta- 
tion and contract the other. The net result is that the small 
scale loops are predominantly of a single polarisation at late 
times. This corresponds to a net current in the direction 
of the cosmic-ray streaming when averaged over the area 
enclosed by ro, i.e. (V x B ■ x) 7^ 0. As the total current en- 
closed by ro increases, the magnetic field Ben on this scale 
must likewise increase. Unlike the small scale fields, however, 
the growth will be independent of orientation. 

To quantify the above simple picture, we combine Q, 
|7| and |8j, together with the equation dA/dt — u x B, 
to give the following expression for the evolution of the fila- 
mentation 



""jcr + ([(«■ V)«]- V)j cr 



(12) 



dt 2 p J v ' dt 

The second two terms on the right hand side of equa- 



tion (121 represent the advection of An with the flow. On 
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small scales, where the velocity gradients are steep, these 
terms dominate over the first term. However, since the con- 
tinued expansion of small scale loops is eventually inhibited 
by neighbouring cavities, (Bell 2004 Reville et al.|2008| , on 
sufficiently large length scales the first term will dominate. 
The ordering of these terms will be verified in section [4] 
Neglecting the last two terms in equation ( 12 1, we find 



the following growth rate for the filamentation instability 



r fi , = 



XB 
P 



"(t) 



t/cr 
P U sh 



1/2 



eBT s 



(13) 



where U C i is the cosmic-ray energy density, 7 m in = Pmin/mc 
the Lorentz factor of the lowest energy cosmic rays driv- 
ing the instability (i.e. those satisfying pmi n c 3> eit s h-A||), 
and r) is a numerical factor that depends on the shape of 
the cosmic-ray spectrum. For a spectrum / oc p~ 4 in the 
momentum interval (p max > p > Pmin), this parameter is 
i] = l/^/ln(p max /p m i n ) . The growth rate is scale indepen- 
dent, and depends only on the root mean square of the 
perpendicular magnetic field enclosed on that scale, as ex- 
pected from the qualitative description above. Since the non- 



resonant mode discussed in Bell (20041 has a growth rate 



that decreases monotonically with increasing wavelength, 
the filamentation must dominate the amplification of mag- 
netic field on some scale. However, for the growth rate of the 
filamentation instability to be sufficiently rapid to influence 
the scattering of high-energy cosmic rays, the mean squared 
magnetic field on small scales must be amplified to values 
well in excess of the ambient field. Thus, the filamentation 
instability can be considered as a bootstrap to the non res- 
onant instability described in |Bell| ( [2004| and |Bell| ( [2005 1. 
The necessary conditions for the filamentation to play an 
important role are discussed in detail in section [4] 

The transfer of magnetic energy from small scales to 
longer wavelengths in the context of diffusive shock accel- 
eration has previously been suggested to occur via an in- 
verse cascade (e.g. |Pelletier et al.|2 006 ; Diamond & Malkov 
20071. Whether this cascade can bridge the large separation 



of scales remains uncertain. The mechanism described here 
presents a different approach where the coupling of the scales 
is mediated by the filamentation. The coupling of small and 
large scale magnetic fields has also been found in simula- 



tions of sheared flows with small scale turbulence (Yousef 



et al. 2008 I in the context of mean-field dynamo theory. 



This approach has also recently been applied to case of pre- 



cursors with an external cosmic-ray current (Bykov et al. 
20TT||Schure fc Bell|20fL . 



3 NUMERICAL SIMULATIONS 

Numerical simulations are performed to verify the analy- 
sis of the previous section. To investigate these processes, 
it is necessary to have a kinetic description of the cosmic- 
rays. A code has been developed similar to that described by 



Zachary fc Cohen ( 1986 1 and Lucek fc Bell ( 2000 ) , where the 



background plasma is treated as an MHD fluid and the cos- 
mic rays are treated using a particle-in-cell (PIC) approach. 
This method is appropriate for modelling plasmas in which 
there exists a large separation between the relevant length 
and time scales associated with the thermal and non-thermal 
particles' kinetics. Unlike full PIC simulations, which solve 



Maxwell's equations directly, our simulations use the mag- 
netic and electric fields determined from the MHD equa- 
tions, i.e.. the electric field is determined from ideal MHD 
E — —u x B, and the displacement current is neglected. 
The cosmic-rays are evolved in these fields by integrating 
the relativistic equations of motion using the Boris method 



(e.g. Birdsall & Langdon 1985 1 and the resulting current 



and charge densities are deposited on the grid at each parti- 



cle update using the approach of Umeda et al. (2003J). Since 
the particles typically evolve on a shorter time-scale than 
the background MHD, it is necessary to integrate particle 
trajectories over several sub-cycles within each MHD update 
( Zachary fc Cohen|1986 l. While this multiplies the effective 
number of particles in the simulation, it can also smooth 
out phase correlations associated with the cosmic-ray cur- 
rent ( |Zachary et al.|1989| [Tucek fc Bell|2000| >. The coupled 
MHD - cosmic-ray equations are solved using a finite- volume 



Godunov scheme, as described in Reville et al. (20081, with 



the self-consistent inclusion of temporal and spatially vary- 
ing current and charge densities j C r, Qcr- 

The simulations are run using a 2D slab symmetric ge- 
ometry, as in the analysis of the previous section, with a 
periodic grid in y— z plane, and the cosmic ray anisotropy di- 
rected out of the simulation plane in the positive ^-direction. 
The particles are initialised as a mono-energetic distribution 
with a net drift velocity. Maintaining gyrotropy in the sim- 
ulations and thus minimising j± due to particle noise, re- 
quires that a large number of particles per cell be used. For 
the results shown in this paper, the particles are initialised 
with 1024 per cell. To minimise the noise in the particle 
distribution at t — 0, the discrete particle momentum vec- 
tor in spherical coordinates, <f>i and fii, are chosen in an 
ordered manner, such that they satisfy the required distri- 
bution globally to a high degree of accuracy. Here p: — +1 
corresponds to the positive x direction. This is achieved 
by choosing the azimuthal angles with uniform spacing as 
<f>i = 2TvM(~iyi/(N - 1) for i = 0.JV, where TV is the 
total number of particles, and M is the number of com- 
plete rotations through 2n required. We found best results 
for M = 3. Having too many or too few revolutions in </> 
at pitch angles /i close to zero can result in a numerically 
introduced anisotropy in the perpendicular direction after 
only a relatively small number of time steps. The particle 
cosines \jh are chosen between fj,i = — 1..1 with decreasing 
spacing Afii such that the total distribution has a net drift 
with the required velocity. Finally the particles are scattered 
in the y — z plane using mixed radix bit-reversed fractions 
( Birdsall fc Langdon||1985 1. 

The two-dimensional magnetic field is initialised by tak- 
ing a sum of Fourier modes in Au sampled from a user- 
defined spectrum, having random wave vectors fc in the 
y — z plane (see e.g. |Giacalone fc Jo kipii 1999). For the 
results shown in Figure [2j we used a 2D Fourier spectrum 
An(fe) oc 1/[1 + (kL c ) 3 ], such that perpendicular magnetic 
field will have a power spectrum peaking close to the length 
L c . Several different forms for the power spectrum have been 
used, and the general results are the same. The computa- 
tional grid isa512x512 square mesh with periodic boundary 
conditions in the y and z directions. The grid resolution is 
Ax — 0.5, where dimensionless length units are chosen such 
that mc/eBo = 1 with Bo the mean magnetic field strength 
out of the plane. The Fourier modes were selected with uni- 
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Figure 2. Time dependent growth of the filamentation. y— 2 plots of Ay, n/no, and p at t = 1,2,3,4. Cosmic rays current is in 

the positive x direction (out of the plane). Details of the simulation parameters are found in the text. 



form logarithmic spacing on the interval 4 < k/2n < 256 
with L c = 16, such that most of the magnetic field struc- 
ture is on scales much smaller than the size of the box. 
The total magnetic field is then calculated by taking the 
curl, and projecting onto the grid using central differenc- 
ing B = Box + V x A«x. This guarantees that the field is 
divergence free. 



The thermal background is initialised at rest with uni- 
form density and pressure. The cosmic-ray current and 
charge densities are also initially uniform, to within noise 
levels on the grid. A rather large shock velocity of u s h = 0.5c 
is used to enhance the anisotropy. If the departure from 
isotropy is too small, a much larger number of particles 
are required to accurately model the perturbations to the 
cosmic-ray current. Since the shock velocity only appears in 
the determination of the cosmic ray anisotropy, the use of the 
non-relativistic MHD equations remains valid. The numeri- 
cal parameters are chosen to represent values that might be 
expected in a young supernova remnant: mn cr /po = 10~ 6 , 
va/c — 5 x 10 -4 . The particle momentum is chosen such 
that the gyroradius is resolved in the simulation box, and 



for the results shown in Figure [2] the particle gyroradius 
is r g — 20(B I Bo)" 1 . We note that with this low minimum 
cosmic ray energy and high shock velocity, the total cosmic 
ray energy density is U CI /poU^ h « 10~ 4 , which is consid- 
erably smaller than what is expected in young supernova 
remnants. These parameters have been chosen such that the 
various physical processes can be easily identified in the sim- 
ulations, although could also indicate that the effect should 
exist even in shocks that are accelerating very inefficiently. 
In reality, for young supernova remnants, the shock veloc- 
ity would be an order of magnitude smaller and the typical 
cosmic ray energy several orders of m agnitude larger, where 
U cr /poU^ h can be as large as 30% (e.g. Malkov & O'C Drury 
2001 1. Simulations with a much higher energy density re- 



sulted in extremely rapid evacuation of zero density cavities 
in the background fluid which provided insufficient time to 
study the growth of magnetic field. 
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Figure 3. Plot of \B±\ at beginning (t = 0) and end (t = 8) of 
the simulation. Regions of strong magnetic field wrap around the 
cosmic ray filaments where the gradients in are steep. Cavities 
are found to develop on a range of scales, in some cases with larger 
loops enclosing smaller loops. 



3.1 Results 

The evolution of Ay, n CT , and p are shown in Figure [2] 
The correlation between the vector potential and cosmic ray 
number density, as represented by n cr , can clearly be seen. 
The cosmic ray filamentation also results in the creation 
and expansion of low density cavities in the thermal plasma. 
There also exists a correlation between the cosmic-ray cur- 
rent j\\ and the potential An, although the features are not 
as sharp as with n cr - As shown in the previous section, the 
correlation between Am and n cr is due to fluctuations in the 
isotropic component fo- The current, on the other hand, 
depends on higher order expansions in the anisotropy of 
the cosmic ray distribution, and is much more difficult to 
capture numerically, but is clearly evident at early times. 
As the simulation progresses, the net streaming of cosmic 
rays is gradually reduced. This is always to be expected in 
simulations of this type (e.g. Lucek & Bell 2000), since the 
work done on the background plasma by the cosmic rays is 
SW = —jci ■ E, i.e. in order to extract energy from the cos- 
mic rays, the background plasma must generate an oppos- 




kLj2n 

Figure 4. Fourier power spectrum of B± at t = 0,2,4,6,8, in 
ascending order for kL/2n < 10, where L is the length of the 
sides of the simulation grid. 



ing electric field. Thus the bulk cosmic-ray drift motion will 
be gradually decelerated, causing the correlation between 
j« and Aii to become weaker as the simulation progresses. 
The opposing electric field will be strongest in the walls 
surrounding the cavities, which can account for the anti- 
correlations found on small scales around regions of large 
A\\ at late times. The correlation on large scales is neverthe- 
less evident even at the end of the simulation, despite a 13% 
reduction in the net streaming velocity. In a real system, 
the cosmic-ray streaming is continually fed by the cosmic- 
ray pressure gradient, an effect that can not be captured 
with our numerical approach. 

The growth of tabletop features in the 2D plots of Au is 
also observed. As discussed in section [2] in two dimensions, 
the value of An cannot increase, but regions where An is 
large spread out to form the features seen in the top row 
of Figure [2] It follows from |6| that the same must hold 
for the cosmic ray current. In a three dimensional situation, 
any gradients in the x direction will act as a source term 
in Q, allowing the possibility of amplifying Ay, resulting 
in continued self-focusing of the cosmic rays. This may have 
important implications for the upstream escape of cosmic 
rays from supernova remnants. We discuss this further in 
section [5] 

The magnetic field at the beginning and end of the sim- 
ulation are shown in Figure |3] As found m previous simu- 
lations, the regions of strongest magnetic field are concen- 



trated in walls surrounding the low density cavities (Bell 



2004 2005 Reville et al. 20081. The cavities are produced 



by the — j cr x Bi force expanding low density cavities with 
the magnetic field frozen in to the background plasma. The 
largest loops appear to be on scales of r ~ L/10 where L 
is the size of the simulation box. To investigate the growth 
of magnetic field on larger scales than this, we perform a 
Fourier analysis of the field. The evolution of the power 
spectrum for B± is shown in Figure [4] It can be seen that 
structures on all wavelengths kL/2-w < 20 grow monoton- 
ically with time. Since a lot of this structure is due to 
the expanding cavities, we focus on the growth of modes 
with kL/2n < 10, which appear to grow almost scale in- 
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dependently, as expected for the filamentation instability. 
The average measured growth rate is ~ 0.25rfu, assuming 
u s h = 0.5c in equation (131. As already mentioned, in the 
simulations, the cosmic-ray drift velocity is reduced as the 
simulation progresses which may account for this reduction 
in comparison with the theoretically predicted growth rate, 
which assumed u sn fixed. The resulting effect on the jit - An 
correlation may also contribute. 

The simulations were terminated when the plasma den- 
sity fell below a certain threshold. In the 3D simulations of 
|Bell| ( |2004| ) , the expanding cavities can merge as field lines 
slide over one another, a feature that cannot be reproduced 
in 2D. Future 3D simulations will allow longer time evolu- 
tion, and address the issue of saturation of the large scale 
magnetic field structure. 

Finally, we emphasise that the primary difficulty in the 
numerical identification of the instability, is choosing pa- 
rameters such that the filamentation instability dominates 
on sufficiently small length scales, in comparison with the 
size of the simulation box. Ultimately, due to the memory 
limitations, we are left with a relatively small dynamical 
range for which a Fourier analysis can be carried out, ap- 
proximately one order of magnitude. In future simulations, 
using a different numerical technique, we hope to extend the 
dynamical range by at least an order of magnitude, and fur- 
thermore, to investigate the role of filamentation in three 
dimensions. 



4 APPLICATION TO DIFFUSIVE SHOCK 
ACCELERATION 

Observations of synchrotron X-ray emission in the vicinity 
of the shock suggest the presence of magnetic fields on the 
order of 100/iG or even larger in several supernova rem- 
nants. Recent observations have also identified the presence 
of a precursor in SN1006, where it is suggested that the im- 
plied field strengths are larger than the typical interstellar 
value (Rakowski et al. 2011), providing tentative evidence 
for field amplification due to the presence of cosmic rays. If 
these large fields are generated via cosmic-ray streaming, the 
mechanism must account for amplification from typical seed 
fields Bism ~ 3 — 5/xG. Numerical simulations of cosmic-ray 
current driven instabilities on small scales suggest non-linear 
amplification of the fields to values B™ s ~ 30-Bism ~ 100/jG 



(Bell 2004 Riquelme & Spitkovsky 2009). This amplification 



is believed to occur far upstream where only the highest en- 
ergy particles are interacting with the interstellar medium 
turbulence ( |Zirakashvili fc Ptuskin|2008||Reville et al.|2009 |. 
The characteristic length scale of the fields produced by the 
non-resonant instability is on too small a scale to reduce 
the mean free path Amfp of the highest energy cosmic rays, 
as required to accelerate them beyond the Lagage-Cesarsky 
limit (Lagagc & Ccsarsky 1983). To increase the accelera- 
tion rate of the highest energy cosmic rays, the diffusion 
coefficient n « Amfpc must be significantly reduced below 
its Bohm limit in the pre-amplified field, where the mean 
free path is equal to the gyroradius Amfp = r g . If the small 
scale fields are indeed amplified to the levels inferred from 
observations far upstream of the shock via the non-resonant 
instability, as we show here, the filamentation instability can 



grow on a sufficiently short time scale to have an appreciable 
influence on the diffusion of the highest energy cosmic rays. 

From the expression for the growth rate given in equa- 
tion | |13| , it can be seen that the filamentation instability 
operates most effectively in strongly amplified small scale 
turbulence, and when the cosmic-rays driving the instabil- 
ity have lower minimum energy. However, if the energy of 
the cosmic-rays driving the filamentation instability is too 
small, the particles will be trapped. This condition, as de- 
rived in section [2] is pc ^> &A||lt s h. Using the fiducial values 

'm. the 



from 

filamentation operates provided 



Bell ( 2004 1 (Equation (21)) fc~ 



2 x 



-Emm > eA\\Ush ~ efe m a X -B±?i s h 

max \ / ^sh 



10 1 



ViooAiGy \2 x io 13 my V 107m / s 



eV.(14) 



TeV 7-ray observations of most historical supernova rem- 
nants provide conclusive evidence for the presence of cosmic 
rays, either protons or electrons, that satisfy this condition. 

To investigate which mechanism determines the trans- 
port properties of the highest energy cosmic rays in the pre- 
cursor, we compare the growth rate of the filamentation in- 
stability to that of the streaming instability given in |Bell| 
( |2004 1, which has a growth rate, 



j 'cr B k 



for fe n 



> k > r„ 1 



(15) 



This expression is equivalent to equation (111 on replacing 



B e /r by kB (cf. |Bell|2005 |. For fc < r g ion-cyclotron res- 



onance takes over and the growth rate steepens oc fc. Hence 
the growth rate r nr can be considered an upper limit for all 
fc. 

Since the magnetic field amplified by the non-resonant 
instability is on too small a scale to effectively scatter the 
highest energy cosmic rays driving the growth, i.e. those 
with E < -Bmax, these particles will continue to gyrate about 
the mean field -Bo. On the scale of the gyroradius of these 
particles fc = e_Boc/-E m ax, the ratio of the growth rate of the 



filamentation instability, equation (13 1, to the non-resonant 
instability is 



£fii 

L nr 



c E n 



(16) 



where -E m i n is the corresponding energy dominating the 
cosmic-ray current. We note here that the small scale fields 
are amplified over a distance much less than the scale- 
height of the precursor at the outer extremity of the precur- 
sor « fc(£; max )/it sh , suggesting that -E max /.E mm should not 
greatly exceed unity. Thus, provided the small scale fields 
can be driven to non-linear values, the filamentation insta- 
bility will play the dominant role in generating the fields 
required to scatter the highest energy cosmic rays in super- 
nova remnants. 

The two growth rates are equal when 



fcr e = 



H 2 



(17) 



where r 9 = Emi n /eBoc. This indicates, to order of magni- 
tude, the length scale above which the filamentation dom- 
inates over the non-resonant mode. Substitution of the pa- 
rameters used in the simulations suggests the transition oc- 
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curs at kL/2n ~ 10, in agreement with what was found. In 



addition, comparing the terms in equation ( 12 1, it is readily 



seen that (171 corresponds to the scale on which the first 



term becomes comparable with the other terms. 

To demonstrate the important role played by the fil- 
amentation instability, we calculate the growth rate using 
typical parameters for young supernova remnants. Assum- 
ing the magnetic field is amplified initially on small length 
scales to a level comparable with those inferred from obser- 
vations, the typical time-scale for growth of magnetic field 
on long wavelengths by the filamentation instability can be 
as short as 



50 



n 

0.2 



0.1 

X V WO^G 



1/2 



10 7 m/s 
E„ 



mm 

10 14 eV 



yrs(18) 



The value of E m [ n driving the growth is the largest uncer- 
tainty. At the onset of the filamentation, provided the lower 
energy cosmic rays satisfy the condition ( |14[ ), the growth 
can be extremely rapid. As the magnetic fields evolve, the 
scale of the filaments becomes comparable to the gyroradius 
of the lower energy particles, such that _E m i n should remain 
large. Saturation may occur when the high energy particles 
become trapped on the self-generated large scale fields. This 
will almost certainly affect the diffusion of cosmic rays and 
may even alter the transport properties of particles at differ- 
ent energies, which can influence the shape of the spectrum 
| |Kirk et al.||1996| . Future simulations in three dimensions 
will help to elucidate this process further. 



5 DISCUSSION 

While the filamentation of photon or high energy electron 
beams in laboratory laser plasma experiments is a well stud- 



ied phenomenon (e.g. Craxton & McCrory 19841, its analogy 
with cosmic rays has been largely overlooked. In this paper, 
it has been demonstrated both analytically and confirmed 
with numerical simulations, that the filamentation of cosmic 
rays is an important process that can occur in the precursors 
of supernova remnants shocks where diffusive shock acceler- 
ation is taking place. 

In addition, we have identified a mechanism for am- 
plifying magnetic field on large length scales as a result of 
the filamentation. The process provides a natural mecha- 
nism to couple the rapid growth of magnetic field on small 
scales, as driven by the non-resonant instability ( Bell|2004 1 , 
to length scales comparable to, or larger than, the gyroradius 
of the particles driving this instability, avoiding the need for 
an inverse-cascade. The growth-time for this instability can 
operate on time scales as short as a few years, provided the 
small scale fields are amplified to a sufficient level. The rea- 
son for the short growth time as compared with previous 
calculations of linear dispersion relations, is that the insta- 
bility develops in non-linear small scale magnetic fields with 
SB 2 /Bo 3> 1, unlike streaming instabilities that typically 
consider the growth of weak Alfvenic perturbations in the 
interstellar medium, where SB 2 / B 2 <C 1. 

This has immediate implications for the maximum en- 
ergy to which cosmic rays can be accelerated at supernova 



remnant shock fronts. The amplification of magnetic turbu- 
lence on all scales, significantly beyond the limits of quasilin- 
ear theory, remains the most likely possibility for accelerat- 
ing cosmic rays above the knee (e.g. Bell & Lucek 2001 Kirk 
& Dendy 2001 ). Using MHD simulations with a constant ex- 
ternal cosmic-ray current, Reville et al. (2008) demonstrated 



that the non-resonant self generated turbulence, reduced the 
diffusion coefficient of test particles significantly below the 
corresponding Bohm value in the pre-amplified field. How- 
ever, due to the limited dynamic range of these simulations, 
the diffusion coefficient for particles with gyroradii larger 
than the typical structures in the field, converged to the 
small-angle scattering limit, where t he mean free path grows 



rapidly with energy Am.f.p. oc E (Reville et al 



2008 



Zi- 



rakashvili & Ptuskin 20081. Thus, the generation of large 



scale field structure is essential to achieve sub-Bohm dif- 
fusion at E > 10 eV. The filamentation instability can 
grow extremely rapidly once the magnetic field perturba- 
tions have been driven to non-linear levels, and may help to 
significantly reduce the mean-free paths of these particles. 
Simulations in 3D, with a significantly larger dynamic range 
are required to confirm this. 

The topic of cosmic rays escaping the source has been 
reinvestigated in recent years in light of developments in 



our understanding of magnetic field amplification (e.g. Ohira 
et al.|2010||Caprioh et al.|2010||Drury|2011| |. These models 
typically assume a one-dimensional, or spherically symmet- 
ric model of cosmic-ray diffusion. Filamentation is almost 
certainly important in this situation, since unlike the early 



one-dimensional trapping models (e.g Kulsrud & Zweibel 
|1975| |Kulsrud]|1979[) , the expansion of low density cavities 



and self-focussed filaments can in fact assist in the escape. 
This most likely occurs only far upstream, since closer to the 
shock, the currents in the filaments are susceptible to beam- 
hose type instabilities. Again, a three-dimensional analysis 
is required to investigate this further, and will be addressed 
in future work. 
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